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Abstract 

The Bloch oscillations in the zigzag array of the optical waveguides are considered. The multiple 
scattering formalism (MSF) is used for the numerical simulation of the optical beam which propagates 
within the array. The effect of the second-order coupling which depends on the geometrical parameters 
of the array is investigated. The results obtained within the MSF are compared with the calculation 
based on the phenomenological coupling modes model. The calculations are performed for the waveguides 
fabricated in alkaline earth boro-aluminosilicate glass sample, which are the most promising for the C-band 
(A ~ 1530 - 1565 nm). 
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I. INTRODUCTION 


The rapid development of technologies for production of optical waveguides provides a possi¬ 
bility of optical simulation of quantum phenomena inherent in the condensed-matter physics jl|, 
such as Bloch oscillations M, Zenner tunneling JslJlo| . Anderson localization ful H] . dynamic 
localization 113, joj]. Periodic arrays of optical waveguides are a special case of photonic crystals. 
Such systems are of interest due to the band structure of their optical spectrum that is similar to 
the band structure in the electronic spectrum of ordinary crystals |jl5, 161. Because of this, the 
behavior of light in photonic crystals is analogous to the behavior of electron wave function in a 
periodic potential of ordinary crystals. This fact is of considerable interest for practical applica¬ 
tions because it allows to steer the optical signal effectively. Besides, it is interesting from a purely 
scientific point of view, since the optical counterparts of quantum phenomena can be observed 
directly. 

Usually, the theoretical model of the light propagation in the optical waveguides array takes 
into account only the next neighbors coupling (NNC). This model is applicable for the particular 
case of the plane arrays, since the high order coupling is too weak. But sometimes the second-order 
coupling (SOC) should be taken into account. The most evident example of a system with the 
significant SOC is the zigzag array (see Fig. [[]). For such a system the SOC depends on the angle 
6 . Varying this angle, we can make the SOC comparable or even greater then the NNC. The SOC 
can infh 


optics 


uence the diffraction processes 


18 


17] and forming of the solitons in the case of the nonlinear 


, 19]- 


In this paper we investigate the influence of SOC on the optical Bloch oscillations. The phe¬ 
nomenon of Bloch oscillations (BO) takes place in the array with the monotonic change of the 
refractive index of the waveguides as one passes from one homogeneous waveguide to another. For 
the case of the plane array, where only the NNC should be taken into account, the light beam path 
possesses an oscillatory form. This effect is widely investigated both experimentally and theoret¬ 
ically in 2-6]. But for the zigzag array, where the SOC is considerable, the path of the optical 
beam takes more complex form. This phenomenon, known as the anharmonic Bloch oscillation 
(ABO), was predicted in 20] and experimentally confirmed in 2lJ. 

Numerical simulation of various phenomena in systems of interacting optical waveguides, in¬ 
cluding BO and ABO, is usually based on a coupled modes model. The parameters of this model 
are typically determined experimentally. 
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Figure 1: Zigzag optical waveguide array. 


In our works 


22 


23] we proposed another method of numerical simulation which requires 


no additional data except of the geometrical and optical parameters of the array (radii of the 
waveguides and distances between them, the refractive indices of the waveguides and environment). 
Our method uses the multiple scattering formalism based on the macroscopic electrodynamics 


approach. In papers 


22 


, |23] we presented the general algorithm for calculating the path of the 


optical beam in an array of waveguides. 

In the present work we generalize the results obtained in 


22 


23]. We apply the multiple 


scattering formalism to the ABO calculation in the zigzag array. The obtained picture of ABO is 
compared with the results of calculation with use of the coupling mode model. For this purpose, 
we derive the analytical formula for the path of the optical beam in the zigzag array, using the 
coupling mode model with the SOC taken into account. Besides, we use the MSF to obtain the 
dependence of the coupling constant of the SOC as function of the angle 6 (see Fig. [[]). 

Recently a new method to fabricate low bend loss femtosecond-laser written waveguides was 
developed [24, 25]. The waveguides are fabricated in alkaline earth boro-aluminosilicate glass 
sample, so they are the most effective for the C-band (A ~ 1530 — 1565 nm). In this work we 
consider a sample that can be fabricated by means of this technology. 

This paper is organized as follows. In Sec. HI] we describe the multiple scattering formalism. In 
subsection III Al we present the algorithm for calculating the path of the optical beam in an array of 
waveguides. In subsection III Bl we derive the equation of the coupling modes model from the MSF 
and present the formulae for the coupling constants. In Sec. IHIl we represent the pictures of ABO 
obtained by use of MSF for different 6 . We compare these results with the paths of the optical 
beam described by the coupling modes model. Besides, in Sec. II I II we present the dependence of 
the SOC coupling constant on the angle 6 . 
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II. MULTIPLE SCATTERING FORMALISM 


Let us consider the array of N parallel infinite cylindrical dielectric waveguides. The axes of the 
waveguides are parallel to the z-axis. The distance between the axes of the adjacent waveguides is 
denoted by a. All the waveguides are assumed to possess the same radius R but different refractive 
indices rij, j being the waveguide number. The refractive index of the environment is n°. The 
permeability of the waveguide material and the environment is unity. 

Suppose that a guided mode with a frequency u) is excited within the array. Then, all the 
components of the electromagnetic held are proportional to e~' lujt+ ' l ^ z . 

Let us consider the held of the guided mode inside of the array. The held of the guided mode is 
finite in any point into the waveguide. So, the electromagnetic held inside of the j-th waveguide 
may be represented in the form 


E;(r) = e 


_ -iujt-\-i(3z 


^ e lm ^ j ( C jm N urfmiPj) - djm M urfmipj) ) , 
m= 0,±1... 


Hj(r) = e -™ t+iPz rij y ( Cjm M UjPm ( Pj ) + d jm TSS UjPm ( Pj )), Pj < R. 

m= 0 ,± 1 ... 


( 1 ) 


Here c Oj = rijUJ, and pj, cf)j are the cylindrical coordinates of the vector r — rj, where r j is the 
coordinates of the axis of the j-th waveguide. The vector cylinder harmonics M Wj p m (pj) and 
are dehned as follows 


% [5 Tfi f3 

{pj ) G r —— d rf! ( ^jPj ) Up 2 drn (Pj ) T 6 2 dm. ( Pj ) ? 


X,- 


*iPj 


~ , . mujj , tujj T , . . 

G r g Jm\^j Pj) T &<f> d rn (X^ pj J, 


X-Pj 


X; 


( 2 ) 


( 3 ) 


where Xj = yju'j — /3 2 , d m (xj P j ) is the Bessel function, and the prime means the derivative with 
respect to the argument XjPj. 

Let us turn to the electromagnetic held outside of the array. This held is the sum of the 
contributions of all the waveguides, 


N N 

e M=E E jM- H(r) = E H A) 

3 =1 3 =1 


( 4 ) 


The contribution induced by the j-th waveguide and vanishing at pj oo may be represented in 
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the form 


Ej(r) = e 


_ —iuit+ifiz 


^ ^ 6 ^ (tijm ^ lj' fUm^Pj') f3m(Pj'i\ j 


m 


H, r = e 


_ — iut+i/3z 0 


£• 


n y e 

m 


im(f)j 




(pj ) 3~ bjm ^uj'/ 3m(Pj ) ) 5 Pj ^ 


( 5 ) 


Here a;' = n°co. The vector cylinder harmonics M Uj pm{Pj) and N ujPm(Pj) are 

771/3 

•Ncj / /3m(P t 7') ^ Pj) ^ (-^ Pj ) “1“ ^2 H m (^X Pj)i 


x! 


* Pi 


IK T / \ m(jJ TT / / \ 26 0 TTt / t \ 

/3m\Pj ) ^9 H m yX pj J ~\~ ~ H m yX Pj J, 


* 2 Pj 


k’ 


( 6 ) 


( 7 ) 


where x 7 = \J a/ 2 — (3 2 , and H m (>dpj) is the Hankel function of the hrst kind. Note that, for 
(3 = 0, Eqs. (H]) and ([5]) transform into the corresponding expressions in Ref. 2^|, however 
different notations are used there. 

Below we consider the simplest approximation to these equations, namely, the zero-harmonic 
approximation. This means that in Eqs. ([I]) and (J5J) only the terms with m = 0 are taken into 
account. In this approximation there are two kinds of the guided modes, namely the transverse 
magnetic (TM) and transverse electric (TE) modes. For the TM-mode bjo = dj 0 = 0, and for the 
TE-mode cqo = Cjq = 0. As an example, let us consider the TM-modes. 


Ej(r) — Cj 

H,(r) 

o 

3 

£ 

Pj < R. 

(8) 

E j( r ) = a 3 N^q {pj), 

H,(r) 

= aj n° M^/qo (pj) 

Pj > R. 

( 9 ) 


Here and below, a 3 and Cj stand for a j0 and c j0 , and the factor e~ lU}t+l P z is omitted. 

To derive the equations that determine the partial amplitudes cq and Cj, one should use the 
boundary conditions on the surface of every waveguide. The boundary conditions on the surface 
of the j-th waveguide connect the held Ej(r), H ; (r) inside the j-th waveguide and the held E(r), 


H(r) outside. In general, there are four independent boundary conditions 
TM-modes and only two boundary conditions are required: 


23| . However, for m = 0 


[E(R/) 1 , 


Ej(Rj) 



[H(R,)]q 


H,(R,) 


( 10 ) 


here Rj is the radius-vector of a point on the surface of the j-th waveguide. 

The boundary conditions on the surface of the j-th waveguide can be expressed in the most 
convenient form if the contributions of all the waveguides to the held outside the array are expressed 
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in terms of the same argument pj. For this propose we apply the following relations: 


N w '0o(pz) ~ Uij(uj,(3) ’Nuipo(pj), 

M u ^o {pi) ~ Uij(u, P) M w '/3o {pj), l j, 


(11) 


where Uij(oj,{3 ) = H a (>c'rij), H 0 being the Hankel function of the first kind and r'i 3 being the 
distance between the axes of the j-th and the Z-tli waveguides. The relations (fTTl) follow from the 


Graf theorem (see 25|) in the zero-harmonic approximation. 
Thus, it follows from (J9]) that 


Ej(r) = a j Ui j (u},P)'N u ip 0 (p j ), H/(r) = aj n° U tj (uj, /3) M^^pj). (12) 


Substituting (TT2l) to (Uj), one gets 

E(r) = cij N ujl p 0 (p j ) + ^ai LZy (w, p) N^o (Pj), 

1*3 

H(r) = dj n° M ul p 0 (pj) + ^ a i n° Uij(u ;, P) M^o {pj)- 

1*3 

So, the boundary conditions (flUl) take the form 


a 3 Ho(x'R) + ^ a z Uij(u, P) J 0 (x'R) = Cj J 0 (xjR), 


1*3 


:n°co' 


n°uj' 


i—rH^x'R) + Uij(uj,P) i—y Jq{h'R) = Cji ^ 1 J^XjR). 

x 1?y 


x' 


x,- 


Eqs. (fl4l) lead to the following system of equations: 


^3 

Cj = Cj((jJ , P) dj. 

Here 

. Ejri J' 0 {xjR) Jq(x'R) — e°Xj Jq(x 3 R) Jq(x'R) 

~ £ oJ^KjR) H' 0 {x!R) - Ejx' J' Q {xjR) H 0 {x'R) ’ 

, Q) = £%{# 0 (xhB) J'(xhB) - H'^R) Jo(ViZ)} 

J ’ £°x, Jq^R) Jo(XjR) — EjX 1 R) 


(13) 


(14) 


(15) 

(16) 

(17) 

(18) 
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A. Method for the optical beam trajectory calculation. 


The system of equations (fT5l) describes the guided modes of the array of waveguides. This system 
possesses the nontrivial solution only if the determinant of the matrix of this system vanishes, 

5ji 


det 


Uji(u,l3) 


= 0. 


(19) 


aj(uj,/3) 

This equation allows to obtain the propagation constants f3 n of the guided modes for the given 
frequency u, n being the number of a guided mode. There are N solutions of Eq. (fT9j) . 


N 


Let aj(/3 n ) be the normalized solution of Eq. (fT51) . ^ \dj(f3 n )\ 2 = 1. The guided mode of the 

3 = 1 

frequency uj is a superposition of the modes with different j3 n \ 


N 


E(f, r) = e-“ Y. C. «^* E M/UrWft), 

n j =1 

N 

H(i, r) = n° e"*" E 0. e‘»"* E “AMMsmM- 


( 20 ) 


3 = 1 


The coefficients C n determine the superposition. 
Let us introduce the modal amplitude 


a j {z) = Y,C n e if}nZ a j (l3 n ). ( 21 ) 

n 

The functions N aJ / / g ri 0 (p J ), M u'p n o(Pj) vanish rapidly as pj increases. So, the field near the j-th 
waveguide is mainly determined by the partial amplitudes cij(/3 n ). So, the modal amplitude a,j(z) 
represents the behavior of the guided modes properly. The coefficients C n are obtained from the 
boundary condition at z = 0: 

*52 C naj (p n ) = aj ( 0 ). ( 22 ) 

n 

The system of equations (1221) allows to obtain the coefficients C n for given Oj(0). 

Below we suppose that the boundary conditions cq (0) possess the form of the Gaussian beam, 

f 7 -_ ?n j2 

aq(0) = e- J~ +ik ° a i. (23) 


This means that the external source approximately illuminates the ends of the waveguides with 
the numbers jo — & < j < jo + ° and the phase difference between the amplitudes taken at the 
ends of the nearest waveguides is k 0 a. 

Thus, the guided mode can be found as follows: 
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1) Calculate numerically the set of propagating constants f3 n using Eq. (115]) ; 

2) Obtain the amplitudes %(/3 n ) for every j3 n using Eq. (fT5j) ; 

3) For the given boundary conditions find the coefficients C n using Eq. 

4) Calculate the function a,j(z) by means of Eq. (1511) . 


B. Derivation of the coupled modes equation 


Suppose that the refractive index is rij = n 0 + Sn x j, Sn -C n 0 . Then, the solutions /3 -°' ) of the 


equation 


dj(uj, f3 


(oh 


= 0 


(24) 


may be represented in the following form: 

Pj 0) = +aj, a</3^ 0) . 
For the case of the weak coupling, (^/3 n — P^ j / P^ 

1 d 1 


(25) 




( 0 ) 


•C 1. Then, one has 

x (a, - Pf). 


aj(Pn) dPcLj(P) 

Then, Eq. (fT5l) leads to the following result: 

(p n - Pj 0) ^j cijiPk) - X] T'y = 0 

where 


(26) 


¥1 


7 lj{Pn) ~ 


Uij (^5 Pn) 


d 1 


dP CLj(ft) 

/3=/3) 0) 


(27) 


(28) 


For the zigzag array we can take into account the first-order and second-order coupling only, 
i.e. I — j ± 1 and l = j ± 2. Due to the symmetry, jjj-i = Jj,j+i and jj,j -2 = lj,j+ 2 - Besides, 
since the waveguides differ insignificantly, one can neglect the dependence of constants and 

lj,j ±2 on j. Then, for the case of weak coupling and insignificant difference of the waveguides one 
can neglect the dependence of 7 ij(P n ) on p n . Thus, Eq. (1271) takes the form 


Pn ~ Pj°j «j(Ai) - 7i pj-i(Pn) + - 72 (a,_ 2 (/3 n ) + a j+ 2 (P n )j = 0. 


(29) 


Here 71 and 72 stand for 7 y,y± 1 and 7 jj ±2 correspondingly. 











Now we can apply the obtained Eq. (1291) to formulate the equation for the function cij(z). It 
follows from (1291) that the function a,j(z) defined by Eq. (12TT) with any set of coefficients C n satisfies 
the equation of the coupled modes. 

a A z ) + 7i ( a i-i( 2 ) + a i+i(-)) + 72 (a>j- 2 {z) + a j+2 (z )) = 0. (30) 


III. PATH OF THE OPTICAL BEAM 


We apply the developed technique for calculating the optical beam in an array represented in 
Fig. [H The parameters taken for the calculation correspond approximately to the parameters of 
the arrays of waveguides reported in work 1261 . 

The wavelength of the laser source in |26| is A = 1550 nm. We take the waveguide radius 
R = 5A = 7750 nm. The refractive index of the environment is n° = 1.4877, and the refractive 
index of the waveguide j — 0 (the center of the array) is n 0 = n° + 5 x 10~ 3 . These parameters also 


approximately corresponds to the experiments reported in 


26]. For our calculations we assume 


that the separation between the waveguides is a = 3 R, and a variation of refractive indices between 
the nearest waveguides is 5n = rij — rij-\ = 5 x 10 -6 . We take the boundary conditions in the 
form ([23]) with a = 4 and ko = 0. 

We produce the calculation for different values of the angle 9 = 50°, 60°, 70°, 180° (the last 
case coincides with the simple plane array). The results of the calculations are represented in Figs. 

H-El 


Every of these figures contains a dotted curve representing the trajectory calculated by means 
of the coupled modes model. The parameter a is found with the help of Eq. (1241) . The coupling 
constants 71 and 72 for all the values 9 are calculated by the use of Eq. (1281) . The values a and 71 do 
not depend on 9 and they are the same for all the four Figs. E1~E1 ol = 16.48 m _1 , 7 x = —58.44 m -1 . 
But 72 depends on the distance b between the second neighbors, so 72 decreases with increasing 9: 
7 2 (50°) = -196.63 m -1 , 7 2 (60°) = -58.44 m" 1 , 7 2 (70°) = -18.64 nr’ 1 , 7 2 (180°) = -0.0277 m” 1 . 
The dependence of the relation 72/71 on the angle 9 is depicted in Fig. El 

Note that for 9 = 60° one has 72 = 71 . This is quite expected result, since for this case b = a. 
Finely, note that for 9 = 180° one has 72/71 ~ 10~ 4 . Thus, we confirm that for the plane array of 
waveguides one can neglect the SOC. 
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Figure 2: 0 = 50° 


IV. CONCLUSION 

In this paper we investigated the anharmonic Bloch oscillations in the zigzag array of the optical 
waveguides. For this purpose we applied the muliple scattering formalism based on the macroscopic 
electrodynamics approach. For the simplicity, we used the zero-harmonic approximation taking 
into account the TM-modes only. On the basis of the MSF, we developed the numerical algorithm 
for calculating the optical beam path for the specified boundary conditions at z = 0. We chose 
the boundary conditions possessing the form of the Gaussian wave packet. It was shown that for 
this case the optical beam path takes the complicated periodic form. 

This result is compared with the optical beam path predicted by means of the coupled modes 
model with the second-order coupling taken into account. This model contains some parameters, 
namely, namely, the difference between the propagation constants of the adjacent waveguides 
a = j3j and the coupling constants 71 and 72 for the first-order and the second-order 
coupling. These parameters were calculated by means of MSF. We showed that the results of 
calculations by two different methods are exactly the same. 
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Figure 3: 0 = 60° 
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Figure 4: 0 = 70° 


11 








70 


s 

o 

nJ 60 


<u 

o 

c 

03 


50 


CA> 


-5 

c 

_o 

os 

bD 

a 

Oh 

o 

S-H 


Oh 


40 

30 

20 

10 


0 



-30 -20 -10 0 10 20 30 


Waveguide 


Figure 5: 8 = 180° 

Besides, we investigated the dependence of the coupling constant of the second-order coupling 
72 on the angle 9. It was shown that 72 decreases quickly as the angle 9 increases. For 9 = 180° the 
value 72 becomes negligible. This result confirms that for the plane array it is possible to produce 
the calculations taking into account the next neighbors coupling only. 
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